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We apply the biorthonormal transfer- matrix renormahzation group (BTMRG) [Phys. Rev. E 83, 
036702 (2011)] to study low-temperature properties of quantum spin chains. Simulation on isotropic 
Heisenberg spin- 1/2 chain demonstrates that the BTMRG outperforms the conventional transfer- 
matrix renormahzation group (TMRG) by successfully accessing far lower temperature unreachable 
by conventional TMRG, while retaining the same level of accuracy. The power of the method is 
further illustrated by the calculation of the low-temperature specific heat for a frustrated spin chain. 
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Quasi-one-dimensional (QID) quantum spin systems 
have been the focus of intensive research for the past 
decades. Quantum fluctuation plays an important role in 
these systems, and powerful non-perturbative analytical 
methods are available (See, for example, Ref. [1]). Re- 
cently, there is a resurgence of interests in QID magnetic 
materials with spin spiral states such as Rb2Cu2Mo30i2 
[2], LiCuV04 0, Li2ZrCu04[4], and LiCu2 02l5|, due to 
their close association with multiferroicity [6| . Typically 
these systems are frustrated and Jordan- Wigner transfor- 
mation approaches can only be applied in limited cases. 
Consequently numerical methods become the major tools 
for understanding these systems. It is well known that 
the density matrix renormahzation group (DMRG) [7| 
is the most powerful numerical method to study ground 
state properties of QID strongly correlated lattice mod- 
els with extremely high precision. DMRG is further de- 
veloped into the transfer matrix renormahzation group 
(TMRG) to study thermodynamics at finite temperature 
by mapping a ID quantum system onto a 2D classical 
counterpart and representing the partition function as 
the trace of powers of quantum transfer matrices (QTMs) 
JSJI . TMRG has been applied to study a variety of quan- 
tum spin chain systems ^, 1iq|. In particular, studies on 
frustrated QID spin chains have predicted several exotic 
quantum phases and very rich phase diagrams. However, 
the behavior at very low temperature remains difficult to 
study [10|. This is because the conventional TMRG suf- 
fers from numerical instabilities at low temperatures due 
to the difficulties in accurately determining the eigenval- 
ues and eigenvectors of the non-Hermitian reduced den- 
sity matrix [8]. Improved TMRG schemes which can ac- 
cess low temperature regime are hence called for. 

In this Letter^ we apply the biorthonormal TMRG 
(BTMRG) method [ll| to accurately determine thermo- 
dynamic quantities of ID spin chains at low tempera- 
tures far below TMRG can reach, while retaining the 



same level of accuracy. BTMRG is built upon a series 
of dual biorthonormal bases for the left and right domi- 
nant eigenvectors of the non-Hermitian QTM [11]. Here, 
the dual biorthonormal bases indicate any two sets of 
vectors {\(y)} and {|/3)} satisfying (a| /3) = 5cx(3- From 
the results of the ID spin-1/2 Heisenberg and frustrated 
J1-J2 models, we believe that BTMRG can reach temper- 
atures that is far below TMRG can reach without losing 
accuracy or suffering from numerical instability. This 
opens up new possibilities to study interesting physics 
such as finite-temperature behaviors near a quantum crit- 
ical point, which are inaccessible previously by other nu- 
merical methods. 

Let us consider the Hamiltonian ilf of a ID quantum 
system of N (even) sites, 
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Using Trotter- Suzuki decomposition [l2| and inserting 
2M complete sets of states {|crp} with site index i and 
Trotter index /c, the partition function can be written as 
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where e = 1/MT and the periodic boundary con- 
ditions along both spatial and Trotter directions 
are assumed. This maps a ID quantum sys- 

tem onto a 2D classical transfer-matrix tensor net- 
work. We define the QTM with length 2M as 
Tm = (^1,2^3,4 •••^2M-i,2m) (^2,3^4,5 •••^2M,i), where 
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FIG. 1. (Color online) (a) A single vertex in the QTM. (b) 
Configuration of the augmented system (environment) block 
Tm (Om), where {\a)} and {|/3)} ({|0} and {|0}) represent 
the reduced dual biorthonormal bases of the current system 
(environment) block. The superblock Tm is formed by con- 
tracting indices s, and t. (c) The enlarged augmented system 
(environment) block Fm+i (Bm+i) is obtained by adding one 
vertex and mapping it onto the new dual biorthonormal bases 
associated with {\a)} and {|/^)}, obtained through renormal- 
ization of Tm • 
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FIG. 2. (Color online) Main steps of RG iteration in BTMRG. 
(a) Computation of the maximal eigenvalue and left and right 
eigenvectors l^^), \(j)) of Tm- (b) The reduced density matrix p^ 
is formed, using \'ij;) and 10) and tracing out the environment. 
Performing SVD on p^ and (c) keeping two reduced sets of 
basis vectors {|a)} and {\b)} from singular vectors of p^ in L 
and R. (d) The bases are biorthonormalized to form a new 
reduced dual biorthonormal bases \{ct)} and {|/3)} for the new 
system block as in Fig. 1(c). 



the site index can be suppressed due to the transla- 
tional invariance. Note that Tm is real-valued but non- 
Hermitian. The maximum eigenvalue and the corre- 
sponding eigenvectors of Tm determine all the physical 
properties in the thermodynamic limit. In practice the 
imaginary-time step s is usually kept fixed. For low tem- 
peratures (M large) the size of Tm is beyond the reach 
of exact diagonalization and DMRG algorithm is applied 
to approximately determine the maximal eigenvalue and 
eigenvectors of Tm [8] . 

The heart of the BTMRG is the bi-orthonormalization 
procedure for any two sets of arbitrary vectors {|a)} 
and {|6)}, where a, 6 = !,••• ,m. Let A = [\a)] and 
B = [\b)] be matrices with the vectors \a) and \b) as 
their columns. Performing a singular value decompo- 
sition (SVD) upon A^B to obtain A^ B = UAV^, we 
can readily obtain a^ dual set of biorthonormal bases 
A = [\a)] = AP and B = [\f3)] = BQ, where P = UA'^^^ 
and Q = ]/A~^/^ represent non-unitary basis transfor- 
mations. In the spirit of DMRG, as shown in Fig. [IJb), 
Tm is partitioned into a system block and an environ- 
ment block (enclosed in dashed frames) together with two 
additional time slices (labeled by a^, a^ and cf^^ct'^). 
At any step, we keep a dual set of biorthonormal bases 
for current system (environment) block labeled by {\a)} 
and 11/3)} ({10} and { I C)}) for the left and right dom- 
inant eigenvectors. Matrix elements of the augmented 
system block Vm are obtained by projecting Tm onto 
the dual biorthonormal bases associated with {|a)} and 
11/3)} as TM{sa(jg,(T'^l3t) = (saa^l Fm Icr'^;/^^), where 
\cr'^j3t) = |cr'£;)(8)|/3)0|t), and similarly for the augmented 
environment block Om onto the bases associated with 



{\i)} and {|C)}- III the case where the Hamiltonian has 
the reflection-symmetry /ii,i+i = /i-i+i,*, Om is simply the 
reflection of Fm, OM(cr£;^t, ^Ccr'^) = FmI^Cct^, cr^;^^), 
such that only the augmented system block Fm needs to 
be calculated and stored. The superblock Tm is built by 
connecting Fm and Om through contracting the indices 
s and t: 



= ^TM{sa(7s,(j'Ep^)^M{crE(.t, seer's). 



(3) 



Figure [2] sketches the main steps of the renormaliza- 
tion group (RG) iteration in BTMRG. We first calculate 
the maximum eigenvalue and the right eigenvectors | (j)) of 
Tm (Fig. [2]^a)). The left eigenvector \il)) can be directly 
read out from |0): i^^^^^^as = ^pcr^Ccrs M- ^^ explicit 
construction of Tm is required since only matrix- vector 
multiplication is involved in determining the maximum 
eigenvalue and eigenvectors, leading to a dramatic reduc- 
tion of memory usage and computation time. The re- 
duced biorthonormal bases {|^)} and {|/3)} for the new 
system block are determined from a non-Hermitian re- 
duced density matrix pg = Tr^j^^ |(/)) ('^l (Fig. Efb)). By 
carrying out SVD upon pg = RT^L^ ^ we keep two reduced 
sets of vectors {\a)} and {\b)} from singular vectors in 
L and R corresponding to the m largest singular val- 
ues in E (Fig. [2fc)). Finally, we perform the prescribed 
biorthonormalization procedure to obtain the dual bases 
{\a)} and {|/3)}. This process is equivalent to applying 
non-unitary basis transformations P and Q upon {\a)} 
and {\b)} (Fig. [2];d)). This completes a cycle of the RG 
steps in BTMRG. The augmented system block is en- 
larged by adding one vertex to Fm, and the enlarged 



Hilbert space is truncated by mapping back to the new 
reduced biorthonormal bases associated with {\a)} and 
{\(3} (Fig.lIJc)). The next cycle is repeated with the new 
system block until the desired M is reached. 

Instead of performing SVD, conventional TMRG diag- 
onalizes pg to obtain two reduced sets of left and right 
eigenvectors, which automatically satisfy the biorthonor- 
mal condition. The BTMRG makes clear this insight 
and generalizes the conventional TMRG to a broad class 
of biorthonormal bases. We note that the choices of the 
biorthonormal bases are not unique [11] and the the steps 
proposed above are crucial in order to exploit symme- 
tries of the Hamiltonian. The advantage of BTMRG over 
the conventional TMRG is now clear. In TMRG, the 
complete diagonalization of the non-Hermitian reduced 
density matrix will inevitably encounter the problem of 
complex eigenvalues and eigenvectors, and introduces nu- 
merical instability that prevents one from reaching very 
low temperatures. In BTMRG, on the other hand, only 
SVD is involved and we obtain real singular values and 
singular vectors. Note that the small singular values ob- 
tained during the biorthonormaliztion procedure can also 
induce instability to the BTMRG algorithm. But this 
can be remedied by replacing the dual corresponding; ba- 
sis vectors associated with the small singular value [11| . 

We demonstrate the power of BTMRG using isotropic 
Heisenberg spin- 1/2 chain whose Hamiltonian reads: 
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The local Hamiltonian /ii,i+i in Eq. (H]) conserves total 
spin thus Qi = Xl/i.(— l)*^^cr^ can be regarded as a good 
quantum number of the QTM [8|. Consequently both 
Tm and the reduced density matrix pg are block-diagonal 
and the SVD and the biorthonormalization can both be 
carried out independently for each subblock for different 
Qi. In addition, the maximum eigenvalue occurs in the 
subblock labeled by gi = |l3(, which results in further 
simplification. The Helmholtz free energy in the ther- 
modynamic limit can be calculated from the maximum 
eigenvalue Aq of Tm as 
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Other thermodynamics quantities can either be calcu- 
lated from the numerical derivatives of the free energy or 
directly using the dominant eigenvectors and eigenvalues 
of the QTM. Systematic errors come from two sources: 
the finiteness of the imaginary-time step £, and trun- 
cation of the reduced basis set. In principle, the first 
source of error has a 0(6:^) correction in the partition 
function. In practice we find that e = 0.05 is adequate 
for the models studied in this work. On the other hand, 
the truncation error can be estimated by the discarded 
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FIG. 3. (Color online) (a) Deviations of the BTMRG free 
energy from the exact results for the isotropic Heisenberg 
model down to temperature T/J = 0.0002 for both m = 80 
and 120. (b) BTMRG and exact results for the specific heat 
Cv. Inset: Low temperature regime where the exact results 
is Cv = (2/3)T. The TMRG arrows indicate the typical tem- 
perature reachable by conventional TMRG. 



weight Wd = 1 — ^^i Ai, where A^'s are the singular val- 
ues in the SVD of the reduced density matrix p^ and m 
is the number of basis states kept by BTMRG. In this 
work, we fixed the discarded weight to be Wd = 10~^^ 
and set up a maximum number of basis states m which 
leads to a fast computation in the first several hundred 
iterations of the RG steps. 

In Fig. [3Ka) we plot the deviation of the BTMRG re- 
sults with m = 80 and 120 from the solution based on the 
Bethe ansatz [l4,ll5|. We find that the BTMRG is very 
stable and can reach temperature down to T/J = 0.0002 
for both 771 's. For a similar tti, conventional TMRG can 
only reach down to T/J = 0.01 due to the numerical 
instability. Furthermore, we note that T/J = 0.0002 is 
not yet the limit of our algorithm, and even lower tem- 
perature can be reached by continuing the RG iteration 
before eventually the numerical instability sets in. More 
importantly, the accuracy is competitive with the con- 
ventional TMRG down to temperature T/J = 0.01 0, 
and the absolute error remains within the order of 10~^ 
down to T/J = 0.0002. At low temperatures the trunca- 
tion error dominates and higher accuracy can be reached 
by keeping more states, which is clearly shown in Fig. [H 
At high temperatures, the Trotter error due to the finite 
imaginary time step dominates. Due to the small length 
of the QTM, the number of states kept by BTMRG to 
keep the discarded weight Wd = 10~^^ may not reach the 
cutoff m. Consequently, the accuracy is less sensitive to 
m. Since the calculation is not variational, the numerical 
results might be larger or smaller than exact values, and 
an zero error crossing will give rise to the cusps in the 
error curve. In Fig. [31(b) we show the results of the spe- 
cific heat C^, and an excellent agreement with the exact 
result is obtained. In the inset, we zoom into the low tem- 
perature regime where the exact result is Cy = (2/3)T 
[8] . The BTMRG results start to deviate from the exact 
value at a much lower temperature compared to conven- 
tional TMRG. The deviation is due to accumulation of 
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T/IJ,I 
FIG. 4. (Color online) Specific heat Cv as a function of 
the temperature for the J1-J2 model down to temperature 
T/\Ji\ — 0.003. Inset: The inflection of the Cv near below 
T /\Ji\ = 0.01 is a result of accumulation of truncation error. 



truncation error and not the numerical instability. 

To examine the scalability of BTMRG, we also study 
the ID frustrated J1-J2 model with the Hamitonian: 
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H — y JiSi • Si-^i -\- J2Si • Si-^2' 



(6) 



We consider the case of ferromagnetic nearest-neighbor 
(NN) interaction (Ji = —1) and antiferromagnetic next- 
nearest-neighbor (NNN) interaction (J2 = 0.4). We 
block two adjacent spins as a superspin to obtain a new 
local Hamiltonian with only NN interaction. Moreover, 
the QTM of this model possesses neither reflection sym- 
metry nor conserved quantum numbers, which makes 
the computation much more numerically demanding. In 
Fig.|4]we calculate the specific heat Cy down to tempera- 
ture T/l Ji| = 0.003 for m = 80 and 120. We find a sharp 
low-temperature anomaly and a broad shoulder at higher 
temperatures which are in good agreement with the re- 
sults of conventional TMRG 0. In the inset of Fig. H 
the low-temperature behavior is magnified. We observe 
an inflection of the Cy curve at very low- temperature. 
For larger m the Cy starts to inflect at lower tempera- 
ture, indicating the inflection is due to the accumulation 
of truncation error and not the numerical stability. By 
increasing the value of m, the BTMRG can still be ac- 
curately pushed to far lower temperature unreachable by 
TMRG without suffering from numerical instability. 

Our BTMRG computations were performed by using 
Matlab on a laptop with an Intel Core i5@ 2.3GHz CPU 
and 4G RAM. It takes about five days to generate a sup- 
perblock with size 2M = 2 • 10^ for m = 120. Since 
BTMRG does not involve complicated eigenvalue de- 
composition of the non-Hermitian p^ as in conventional 
TMRG, no special numerical tricks are necessary, and 
the computational complexity is significantly reduced. 

In summary, we apply BTMRG method to accurately 
determine the thermodynamic properties of ID quantum 
spin chains at low temperatures. Due to the improved 
numerical stability of the BTMRG, the simulations can 



access temperatures several orders of magnitude lower 
than the conventional TMRG, while retaining the same 
level of accuracy. This work suggests that the BTMRG is 
a promising numerical method to investigate the low tem- 
perature thermodynamics of QID quantum systems. For 
example, it may become possible to study the behavior 
near a quantum critical point, the transition from high- 
temperature regime into the quantum critical regime, and 
make direct comparison with experiments. In addition it 
may also be possible to extend the long-time limit in the 
studies of real-time dynamics [17| . 

We acknowledge the support by NSC in Taiwan 
through Grants No. 100-2115-M-232-001 (Y.K. Huang.), 
98-2112-M-007-010 (P. Chen), 100-2112-M-002-013-MY3 
(Y.J. Kao), and by NTU Grant numbers 10R80909-4 
(Y.J. Kao). Travel support from NCTS in Taiwan is 
also acknowledged. Y.K. Huang would also like to give 
special thanks to Prof. Tao Xiang for useful discussions. 



ykln@mail.njtc.edu.tw| 



pcchen@phys.nthu.edu.tw| 
'yjkao@phys.ntu.edu.tw| 



[1] T. Giamarchi, ^Quantum Physics in One Dimension ^ 

[2] 



[3] 



[4] 



[5] 

[6] 
[7] 



(Oxford University Press, USA, 2004). 
M. Hase, H. Kuroe, K. Ozawa, 
H. Kitazawa, G. Kido, and 

|Phys. Rev. B 70, 104426 (2004) 



O. 
T. 



Suzuki, 
Sekine, 



M. Enderle, C. Mukherjee, B. Fak, R. K. Kremer, J.-M. 
Broto, H. Rosner, S.-L. Drechsler, J. Richter, J. Malek, 
A. Prokofiev, W. Assmus, S. Pujol, J.-L. Raggazzoni, 
H. Rakoto, M. Rheinstadter, and H. M. R0nnow, 
EPL ( Europhysics Letters) 70, 237 (2005) , 
's.-L. Drechsler, O. Volkova, A. N. Vasiliev, N. Tris- 
tan, J. Richter, M. Schmitt, H. Rosner, J. Malek, 
R. Klingeler, A. A. Zvy agin, and B. Biichner, 
"Phys. Rev. Lett. 98, 07720272007) 
S. Park, Y. J. Choi, C. L. Zhang, 
Phys. Rev. Lett. 98, 057601 (2007), 
S.-W. Cheong and M. Mostovoy, Nat Mater 6, 13 (2007) 
S. R. White, Phys. Rev. Lett. 69, 2863 (1992)^ 
Phys. Rev. B 48, 10 345 (1993); U. Schollwock, 

Rev. Mod. Phys. 77 , 259 (2005) ; 
Advanc es in Physics 55, 477 (2006) 
Xiang, 



and S.-W. Cheong, 



K. A. Hallberg, 



[8] R. J. Bursill, T^ 



and G. A. Gehring, 



Journal of Physics: Condensed Matter 8, L583 (1996)J 



[9] 

[10] 

[11] 
[12] 



)J 
X. Wang an d T. Xiang, |Phys. Rev. B 56, 5061 (1997J1 

T. Xiang, P hys. Rev. B 58, 9142 (1998)] F. Naef, 

X. Wa ng, X. Zotos, and W. von der Linden, 

Phys. Rev. B 60, 359 (1999); J. Sirker and A. Kliimper, 

EPL (Europhysics Letters) 60 , 262 (2002j| J. Sirker, 

Phys. Rev. B 69, 104428 (20047] 

H. T. Lu, Y. J. Wang, S. Qin, and 



Xiang, 
Sota 



Phys. Rev. B 74, 134425 (2006) | 



and 



T. 



Tohy ama, 
Journal of Physi cs: Conf erence Series 200, 012l91 (20l6)] 
Y.-K. Huan^; ^Phys. Rev. E 83, 036702 (2011); 

J. Stat. Mech.: Theory Exp. 2011, P07003 (2011), 
H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1958); 



M. Suzuki, Common. Math. Phys. 51, 183 (1976)] Physics, Vol. 645 (Springer Berhn / Heidelberg, 2004) 

10.1007/BF01609348; „Ph ys. Rev. B 31, 2957 ( 198j)j pp. 349-379. 

l^^] K. Nomura and M. Yamada, [16] J. Sirker, Transfer matrix approach to thermodynam- 



Phys. Rev. B 43, 8217 (1991) ics and dynamics of one- dimensional quantum systems^ 



[14] A. Kliimper, Zeitschrift fiir Physik B Condensed Matter 91, 50 7 (3 a9D), thesis, Universitat Dortmund (2002). 

10.1007/BF01316831. ^ ^^^^ [17] J. Sirker and A. Kliimper, 

[15] A. Kliimper, in Quantum Magnetism , Lecture Notes in Phys. Rev. B 71, 241101 (2005) _ 



